#define SOL_ALL_SAFETIES_ON 1
#include <sol/sol.hpp>
#include "matrix.hpp"
#include "eigen.h"
#include "eigen_spectra.h"

#include <iostream>
#include <list>
#include <cmath>
#include <complex>

typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> EigenMatrix;
typedef Eigen::MatrixXcd EigenCMatrix;
typedef Eigen::VectorXcd EigenVectorXcd;
typedef Eigen::VectorXd EigenVectorXd;

typedef Eigen::Matrix<std::string, Eigen::Dynamic, Eigen::Dynamic> EigenSMatrix;


namespace sol {
        //template <>
        //struct is_automagical<EigenMatrix> : std::false_type {};
        //template <>
        //struct is_automagical<EigenVectorXcd> : std::false_type {};
}

// automatically bound to tostring(obj) [ __tostring ]

namespace Eigen{
	std::string serialize( const EigenSMatrix& mat )
	{
		std::string str = "nl.new_eigen_smatrix(";
		str += std::to_string(mat.rows()) + "," + std::to_string(mat.cols()) + ",{";
		for(int ii=0; ii < mat.rows(); ii++){
			for(int jj=0; jj < mat.cols(); jj++){
				str += "'" + mat(ii,jj) + "';";
			}
		}
		str += "}); ";
		return str;
	}
}

template< class T >
std::ostream& operator<<(std::ostream& os, const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& right) {

	os << "ematrix" <<std::endl;
	if (right.size() == 0) {
		os << "{ ematrix! empty }";
		return os;
	}

	os << "{ ematrix " << right.size() << " | ";
	for(int ii=0; ii <right.size(); ii++){
		os << right(ii);
		if(ii!=right.size()-1)
			os << ", " ;
	}
	os << " }";
	return os;
}

std::ostream& operator<<(std::ostream& os, const EigenVectorXcd& right) {

	os << "eigen_cvector" <<std::endl;
	if (right.size() == 0) {
		os << "{ eigen_cvector! empty }";
		return os;
	}

	os << "{ eigen_cvector " << right.size() << " | ";
	for(int ii=0; ii <right.size(); ii++){
		os << right(ii);
		if(ii!=right.size()-1)
			os << ", " ;
	}
	os << " }";
	return os;
}

std::ostream& operator<<(std::ostream& os, const EigenVectorXd& right) {

	os << "eigen_vector" <<std::endl;
	if (right.size() == 0) {
		os << "{ eigen_vector! empty }";
		return os;
	}

	os << "{ eigen_vector " << right.size() << " | ";
	for(int ii=0; ii <right.size(); ii++){
		os << ", " << right(ii);
	}
	os << " }";
	return os;
}

EigenMatrix operator+(double left,const EigenMatrix& right) {
	EigenMatrix res = right;
	res = right*left;
	return res;
}
EigenMatrix operator-(double left,const EigenMatrix& right) {
	EigenMatrix res = right;
	res = right*left;
	return res;
}
namespace LuaNumeric {


template<class T>
sol::object getter(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& vec, sol::table key, sol::this_state L) {
		// we use stack_object for the arguments because we know
		// the values from Lua will remain on Lua's stack,
		// so long we we don't mess with it
		// String keys failed, check for numbers
	sol::object key1 = key[1];
	sol::object key2 = key[2];

	auto maybe_numeric_key1 = key1.as<sol::optional<int>>();
	auto maybe_numeric_key2 = key2.as<sol::optional<int>>();
		if (maybe_numeric_key1&&maybe_numeric_key2) {
			int n1 = *maybe_numeric_key1;
			int n2 = *maybe_numeric_key2;
			if( n1 <= vec.rows() && n1 >= 1
				&& n2 <= vec.cols() && n2 >= 1)
				return sol::object(L, sol::in_place, vec.coeff(n1-1, n2-1));
		}
		// No valid key: push nil
		// Note that the example above is a bit unnecessary:
		// binding the variables x and y to the usertype
		// would work just as well and we would not
		// need to look it up here,
		// but this is to show it works regardless
		return sol::object(L, sol::in_place, sol::lua_nil);
}
template<class T,class T1>
void setter(Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& vec, sol::table key, T1 val) {
		// we use stack_object for the arguments because we know
		// the values from Lua will remain on Lua's stack,
		// so long we we don't mess with it
		// String keys failed, check for numbers
	sol::object key1 = key[1];
	sol::object key2 = key[2];

	auto maybe_numeric_key1 = key1.as<sol::optional<int>>();
	auto maybe_numeric_key2 = key2.as<sol::optional<int>>();
	if (maybe_numeric_key1&&maybe_numeric_key2) {
		int n1 = *maybe_numeric_key1;
		int n2 = *maybe_numeric_key2;
		if( n1 <= vec.rows() && n1 >= 1
			&& n2 <= vec.cols() && n2 >= 1){
			vec.coeffRef(n1-1, n2-1) = val;
			return;
		}
	}
	throw "Eigen::Matrix out of range or index error";
}

template<class T>
sol::object vgetter(const Eigen::Matrix<T, Eigen::Dynamic, 1>& vec, sol::object key, sol::this_state L) {
		// we use stack_object for the arguments because we know
		// the values from Lua will remain on Lua's stack,
		// so long we we don't mess with it
		// String keys failed, check for numbers

	auto maybe_numeric_key = key.as<sol::optional<int>>();
		if (maybe_numeric_key) {
			int n1 = *maybe_numeric_key;
			if( n1 <= vec.size() && n1 >= 1 )
				return sol::object(L, sol::in_place, vec.coeff(n1-1));
		}
		// No valid key: push nil
		// Note that the example above is a bit unnecessary:
		// binding the variables x and y to the usertype
		// would work just as well and we would not
		// need to look it up here,
		// but this is to show it works regardless
		return sol::object(L, sol::in_place, sol::lua_nil);
}
template<class T>
void vsetter(Eigen::Matrix<T, Eigen::Dynamic, 1>& vec, sol::object key, T val) {
		// we use stack_object for the arguments because we know
		// the values from Lua will remain on Lua's stack,
		// so long we we don't mess with it
		// String keys failed, check for numbers
	auto maybe_numeric_key = key.as<sol::optional<int>>();
		if (maybe_numeric_key) {
			int n1 = *maybe_numeric_key;
			if( n1 <= vec.size() && n1 >= 1 )
				vec.coeffRef(n1-1) = val;
				return;
		}
	throw "Eigen::Matrix out of range or index error at vsetter";
}
	typedef std::complex<double> complex;
	void open_luaEigen(sol::table& module) {
		module.set_function("new_eigen_matrix", sol::overload(
				[]() { return EigenMatrix(); }
				,[](int m, int n) { return EigenMatrix(m,n); }
				)
		);
		module.set_function("new_eigen_cmatrix", sol::overload(
						[]() { return EigenCMatrix(); }
						,[](int m, int n) { return EigenCMatrix(m,n); }
						)
				);
		module.set_function("new_eigen_smatrix", sol::overload(
						[]() { return EigenSMatrix(); }
						,[](int m, int n) { return EigenSMatrix(m,n); }
						,[](int m, int n,sol::table t) {
							EigenSMatrix mat(m,n);
							if( auto res = t.as< sol::optional< std::vector<std::string> > > () )
							{
								if((*res).size() == m*n ){
									for(int ii=0;ii<m; ii++){
										for(int jj=0; jj < n; jj++){
											mat(ii,jj) = (*res)[ii*n+jj];
										}
									}
								}else{
									throw "Data size mismatch";
								}
							}
								return mat;
						}
						)
				);
		module.set_function("new_eigen_vectorxd", sol::overload(
				[]() { return EigenVectorXd(); }
				,[](int m) { return EigenVectorXd(m); }
				)
		);
		module.set_function("new_eigen_vectorxcd", sol::overload(
				[]() { return EigenVectorXcd(); }
				,[](int m) { return EigenVectorXcd(m); }
				)
		);

		module.set_function("deigs", sol::overload(
					[](Eigen::MatrixXd& A,Eigen::MatrixXd& B,int nvec) {
						std::tuple<int,Eigen::VectorXd,Eigen::MatrixXd> res;
						std::get<0>(res) = solve_deig_spectra(A,
													  B,
													  std::get<1>(res),
													  std::get<2>(res),
													  nvec);
						return std::move(res);
					}
				)
		);

		module.set_function("eig", sol::overload(
					[](Eigen::MatrixXd& A,Eigen::MatrixXd& B, const std::string& strSort) {
						std::tuple<int,Eigen::VectorXcd,Eigen::MatrixXcd> res;
						std::get<0>(res) = solve_geig( A,B,
													   std::get<1>(res),
													   std::get<2>(res));
						sortEigenVectorByValues(std::get<1>(res),std::get<2>(res),strSort);
						return std::move(res);
					},
					[](Eigen::MatrixXd& A,Eigen::MatrixXd& B) {
						std::tuple<int,Eigen::VectorXcd,Eigen::MatrixXcd> res;
						std::get<0>(res) = solve_geig( A,B,
													   std::get<1>(res),
													   std::get<2>(res));
						return std::move(res);
					},

					[](Eigen::MatrixXd& A) {
						std::tuple<int,Eigen::VectorXcd,Eigen::MatrixXcd> res;
						std::get<0>(res) = solve_eig( A,
													  std::get<1>(res),
													  std::get<2>(res));
						return std::move(res);
					},
					[](Eigen::MatrixXd& A, const std::string& strSort) {
						std::tuple<int,Eigen::VectorXcd,Eigen::MatrixXcd> res;
						std::get<0>(res) = solve_eig( A,
													  std::get<1>(res),
													  std::get<2>(res));
						sortEigenVectorByValues(std::get<1>(res),std::get<2>(res),strSort);
						return std::move(res);
					},
					[](Eigen::MatrixXcd& A) {
						std::tuple<int,Eigen::VectorXcd,Eigen::MatrixXcd> res;
						std::get<0>(res) = solve_ceig( A,
													  std::get<1>(res),
													  std::get<2>(res));
						return std::move(res);
					},
					[](Eigen::MatrixXcd& A, const std::string& strSort) {
						std::tuple<int,Eigen::VectorXcd,Eigen::MatrixXcd> res;
						std::get<0>(res) = solve_ceig( A,
													  std::get<1>(res),
													  std::get<2>(res));
						sortEigenVectorByValues(std::get<1>(res),std::get<2>(res),strSort);
						return std::move(res);
					}
				)
		);

		module.set_function("eigh", sol::overload(
					[](Eigen::MatrixXd& A,Eigen::MatrixXd& B, sol::object key ) {
						std::tuple<int,Eigen::VectorXcd,Eigen::MatrixXcd> res;
						bool bSort = false;
						auto maybe_bool_key = key.as<sol::optional<bool>>();
						if (maybe_bool_key)
							bSort = *maybe_bool_key;

						std::get<0>(res) = solve_gheig( A,B,
													   std::get<1>(res),
													   std::get<2>(res));

						return std::move(res);
					},
					[](Eigen::MatrixXcd& A,Eigen::MatrixXcd& B) {
						std::tuple<int,Eigen::VectorXcd,Eigen::MatrixXcd> res;
						std::get<0>(res) = solve_gheig( A,B,
													   std::get<1>(res),
													   std::get<2>(res));

						return std::move(res);
					},
					[](Eigen::MatrixXd& A) {
						std::tuple<int,Eigen::VectorXcd,Eigen::MatrixXcd> res;
						std::get<0>(res) = solve_heig( A,
													   std::get<1>(res),
													   std::get<2>(res));

						return std::move(res);
					},
					[](Eigen::MatrixXcd& A) {
						std::tuple<int,Eigen::VectorXcd,Eigen::MatrixXcd> res;
						std::get<0>(res) = solve_heig( A,
													   std::get<1>(res),
													   std::get<2>(res));

						return std::move(res);
					}
				)
		);

		module.new_usertype<EigenVectorXcd>("EigenVectorXcd"
				,sol::no_constructor
				,"resize", (void(EigenVectorXcd::*)(EigenVectorXcd::Index)) &EigenVectorXcd::resize
				,"get", [](EigenVectorXcd* m, int nm )->complex& { return (*m)(nm-1);}
				,"set", sol::overload(
						[](EigenVectorXcd& m, int nm,const complex& v){ m(nm-1) = v;}
						,[](EigenVectorXcd& m, int nm,double v){ m(nm-1) = v;}
				)
				,"set_ones",[](EigenVectorXcd& m){
						m.setOnes();
				}
				,"set_zero",[](EigenVectorXcd& m){
						m.setZero();
				}
				,sol::meta_function::index, &vgetter<complex>
				,sol::meta_function::new_index,  sol::overload(
						[](EigenVectorXcd& m, int nm,const complex& v){ m(nm-1) = v;}
						,[](EigenVectorXcd& m, int nm,double v){ m(nm-1) = v;}
				)
		);
		module.new_usertype<EigenVectorXd>("EigenVectorXd"
				,sol::no_constructor
				,"resize", (void(EigenVectorXd::*)(EigenVectorXd::Index)) &EigenVectorXd::resize
				,"get", [](EigenVectorXd* m, int nm )->double& { return (*m)(nm-1);}
				,"set", [](EigenVectorXd* m, int nm,const double& v){ (*m)(nm-1) = v;}
				,"set_ones",[](EigenVectorXd& m){
						m.setOnes();
				}
				,"set_zero",[](EigenVectorXd& m){
						m.setZero();
				}
				,"size",[](EigenVectorXd& m){
						m.size();
				}
				,"to_std_vector",[](EigenVectorXd& m){
					std::vector<double> vec(m.size());
					for(int ii=0;ii<m.size();ii++){
						vec[ii] = m(ii);
					}
					return vec;
				}
				,sol::meta_function::index, &vgetter<double>
				,sol::meta_function::new_index, &vsetter<double>

		);

		module.new_usertype<EigenMatrix>("EigenMatrix"
				,sol::no_constructor
				//sol::constructors<EigenMatrix(EigenMatrix::Index, EigenMatrix::Index)>(),
				,"resize", (void(EigenMatrix::*)(EigenMatrix::Index, EigenMatrix::Index)) &EigenMatrix::resize
				//,sol::meta_function::addition, sol::resolve< EigenMatrix(double, const EigenMatrix&)>(::operator+)
				//,sol::meta_function::subtraction, sol::resolve< EigenMatrix(double, const EigenMatrix&)>(::operator-)
				,"get", [](EigenMatrix* m, int nm,int nn )->double& { return (*m)(nm-1,nn-1);}
				,"set", [](EigenMatrix* m, int nm,int nn, double v){ (*m)(nm-1,nn-1) = v;}

				,"set_row", [](EigenMatrix& m, int nm,double v){
					if(nm<1 || nm > m.rows()){
						throw " out of range ";
					}
					for( int ii=0; ii < m.cols(); ii++){
						m(nm-1, ii) = v;
					}
				}
				,"set_col", [](EigenMatrix& m, int nm,double v){
					if(nm<1 || nm > m.cols()){
						throw " out of range ";
					}
					for( int ii=0; ii < m.rows(); ii++){
						m(ii,nm-1) = v;
					}
				}

				,sol::meta_function::index, &getter<double>
				,sol::meta_function::new_index, &setter<double,double>
				,sol::meta_function::multiplication, sol::overload(
							[](const EigenMatrix& a, const EigenMatrix& b){
								EigenMatrix c = a*b;
								return std::move(c);
							}

				)
				,sol::meta_function::addition, sol::overload(
							[](const EigenMatrix& a, const EigenMatrix& b){
								EigenMatrix c = a+b;
								return std::move(c);
							}
							//,sol::resolve< EigenMatrix(double, const EigenMatrix&)>(::operator+)
				)
				,sol::meta_function::subtraction, sol::overload(
							[](const EigenMatrix& a, const EigenMatrix& b){
								EigenMatrix c = a-b;
								return std::move(c);
							}
				)
				,"rows",[](EigenMatrix& M){
					return M.rows();
				}
				,"cols",[](EigenMatrix& M){
					return M.cols();
				}
				,"row",[](EigenMatrix& M,int idx,sol::this_state L){
					if(idx<1 || idx > M.rows())
						throw "out of range";
					return EigenVectorXd( M.row(idx-1) );
				}
				,"col",[](EigenMatrix& M,int idx,sol::this_state L){
					if(idx<1 || idx > M.cols())
						throw "out of range";
					return EigenVectorXd( M.col(idx-1) );
				}
				,"inverse",[](EigenMatrix& M){
					return EigenMatrix(M.inverse());
				}
				,"set_identity",[](EigenMatrix& m){
						m.setIdentity();
				}
				,"set_ones",[](EigenMatrix& m){
						m.setOnes();
				}
				,"set_zero",[](EigenMatrix& m){
						m.setZero();
				}

				,"set_random",[](EigenMatrix& m){
						m.setRandom();
				}
				,"solve_lu", sol::overload(
						[](EigenMatrix& M, const Eigen::VectorXd& b, Eigen::VectorXd& x){
							int res = solve_lu(M,b,x);
								return res;
						}
						,[](EigenMatrix& M, const EigenMatrix& B, EigenMatrix& X){
							int res = solve_lu(M,B,X);
							return res;
						}
				)
				,"solve_eig", sol::overload(
						[](EigenMatrix& M, Eigen::VectorXd& v, EigenMatrix& E){
							//int res = solve_eig(M,v,E);
							//return res;
						}
				)
				,sol::meta_function::multiplication, sol::overload(
					[](const EigenMatrix& a, const EigenMatrix& b){
						EigenMatrix c = a*b;
						return std::move(c);
					}

				)
				,sol::meta_function::addition, sol::overload(
							[](const EigenMatrix& a, const EigenMatrix& b){
								EigenMatrix c = a+b;
								return std::move(c);
							}
							//,sol::resolve< EigenMatrix(double, const EigenMatrix&)>(::operator+)
				)
				,sol::meta_function::subtraction, sol::overload(
							[](const EigenMatrix& a, const EigenMatrix& b){
								EigenMatrix c = a-b;
								return std::move(c);
							}
				)
		);


		module.new_usertype<EigenCMatrix>("EigenCMatrix"
				,sol::no_constructor
				,"resize", (void(EigenCMatrix::*)(EigenCMatrix::Index,EigenCMatrix::Index)) &EigenCMatrix::resize
				,"get", [](EigenCMatrix* m, int nm,int nn )->complex& { return (*m)(nm-1,nn-1);}
				,"set", [](EigenCMatrix* m, int nm,int nn, const complex& v){ (*m)(nm-1,nn-1) = v;}
				,sol::meta_function::index, &getter<complex>
				,sol::meta_function::new_index, sol::overload(
						 &setter<std::complex<double>,std::complex<double>>
						,&setter<std::complex<double>,double>
				)

				,"set_row", sol::overload(
					[](EigenCMatrix& m, int nm,double v){
						if(nm<1 || nm > m.rows()){
							throw " out of range ";
						}
						for( int ii=0; ii < m.cols(); ii++){
							m(nm-1, ii) = v;
						}
					},
					[](EigenCMatrix& m, int nm, const std::complex<double>& v){
							if(nm<1 || nm > m.rows()){
								throw " out of range ";
							}
							for( int ii=0; ii < m.cols(); ii++){
								m(nm-1, ii) = v;
							}
					}
				)
				,"set_col",sol::overload(
					[](EigenCMatrix& m, int nm,double v){
						if(nm<1 || nm > m.cols()){
							throw " out of range ";
						}
						for( int ii=0; ii < m.rows(); ii++){
							m(ii,nm-1) = v;
						}
					},
					[](EigenCMatrix& m, int nm,const std::complex<double>& v){
						if(nm<1 || nm > m.cols()){
							throw " out of range ";
						}
						for( int ii=0; ii < m.rows(); ii++){
							m(ii,nm-1) = v;
						}
					}
				)
/**/
				,"rows",[](EigenCMatrix& M){
					return M.rows();
				}
				,"cols",[](EigenCMatrix& M){
					return M.cols();
				}
				,"row",[](EigenCMatrix& M,int idx,sol::this_state L){
					if(idx<1 || idx > M.rows())
						throw "out of range";
					return EigenVectorXcd( M.row(idx-1) );
				}
				,"col",[](EigenCMatrix& M,int idx,sol::this_state L){
					if(idx<1 || idx > M.cols())
						throw "out of range";
					return EigenVectorXcd( M.col(idx-1) );
				}
				,"inverse",[](EigenCMatrix& M){
					return M.inverse();
				}
				,"set_identity",[](EigenCMatrix& m){
						m.setIdentity();
				}
				,"set_ones",[](EigenCMatrix& m){
						m.setOnes();
				}
				,"set_zero",[](EigenCMatrix& m){
						m.setZero();
				}
				,"set_random",[](EigenCMatrix& m){
						m.setRandom();
				}
				,"solve_lu", sol::overload(
						[](EigenCMatrix& M, const Eigen::VectorXcd& b, Eigen::VectorXcd& x){
							int res = solve_lu(M,b,x);
								return res;
						}
						,[](EigenCMatrix& M, const EigenCMatrix& B, EigenCMatrix& X){
							int res = solve_lu(M,B,X);
							return res;
						}
				)
				,"solve_eig", sol::overload(
						[](EigenCMatrix& M, Eigen::VectorXcd& v, EigenCMatrix& E){
							int res = solve_ceig(M,v,E);
							return res;
						}
				)
				,sol::meta_function::multiplication, sol::overload(
							[](const EigenCMatrix& a, const EigenCMatrix& b){
								EigenCMatrix c = a*b;
								return std::move(c);
							}

				)
				,sol::meta_function::addition, sol::overload(
							[](const EigenCMatrix& a, const EigenCMatrix& b){
								EigenCMatrix c = a+b;
								return std::move(c);
							}
							//,sol::resolve< EigenMatrix(double, const EigenMatrix&)>(::operator+)
				)
				,sol::meta_function::subtraction, sol::overload(
							[](const EigenCMatrix& a, const EigenCMatrix& b){
								EigenCMatrix c = a-b;
								return std::move(c);
							}
				)
		);

		module.new_usertype<EigenSMatrix>("EigenSMatrix"
				,sol::no_constructor
				,"resize", (void(EigenSMatrix::*)(EigenSMatrix::Index,EigenSMatrix::Index))
									&EigenSMatrix::resize
				,"get", [](EigenSMatrix* m, int nm,int nn )->std::string& { return (*m)(nm-1,nn-1);}
				,"set", [](EigenSMatrix* m, int nm,int nn, const std::string& v){ (*m)(nm-1,nn-1) = v;}
				,"rows",[](EigenSMatrix& M){
					return M.rows();
				}
				,"cols",[](EigenSMatrix& M){
					return M.cols();
				}
				,sol::meta_function::index, &getter<std::string>
				,sol::meta_function::new_index, &setter<std::string,std::string>
		);

	}
} // namespace my_object




